from general_neuro.conductance_neuron import Conductance_neuron
from general_neuro.alpha_synapse import Alpha_synapse
from scipy.integrate import odeint
from genutils.io import raw_data_columns
from scipy import linspace,arange,linalg,matrix,vectorize,array,hstack,vstack,array
from pylab import plot,show,subplot,figure
import pylab

#----------------------------------------------------------------------------------
# some synaptic input
#----------------------------------------------------------------------------------
holding_V = [-68,-80,-60,-50.0] # in mV
data_columns = [1,3,5,7]

data = raw_data_columns("09-10-08-vclamp-avg.DAT")
data = array(data)

t_imc   = data.T[0][5040:5700]-500.00
t_sgp   = data.T[0][20040:20700]-2000.00
imc = []
sgp = []
for col in data_columns:
    imc.append(data.T[col][5040:5700])
    sgp.append(data.T[col][20040:20700])


def sfc(c1,c2,E_rev1,E_rev2,hp1,hp2):
    # solve for g1 and g2 given voltage clamp info
    a = matrix([[E_rev1-hp1, E_rev2-hp1],[E_rev1-hp2, E_rev2-hp2]])
    b = matrix([[c1],[c2]])
    x = linalg.solve(a,b)
    return tuple(-1*x)

sfc_vectorized = vectorize(sfc)

def calc_conductances(rev_pot):
    
    for i in range(len(pairs)):
        j,k = pairs[i]
        g1[i],g2[i] = sfc_vectorized(imc[j],imc[k],rev_pot,-79.0,holding_V[j],holding_V[k])
        g1_sgp[i],g2_sgp[i] = sfc_vectorized(sgp[j],sgp[k],rev_pot,\
                                -79.0,holding_V[j],holding_V[k])

    # calc vstacks
    for i in range(len(pairs)):
        g1_v = vstack(g1).T
        g2_v = vstack(g2).T
        g1_sgp_v = vstack(g1_sgp).T
        g2_sgp_v = vstack(g2_sgp).T
    # calc averages
    g1_std = []
    g2_std = []
    g1_sgp_std = []
    g2_sgp_std = []
    for i in range(len(g1[0])):
        g1_avg.append(g1_v[i].sum()/len(g1))
        g2_avg.append(g2_v[i].sum()/len(g2))
        g1_sgp_avg.append(g1_sgp_v[i].sum()/len(g1_sgp))
        g2_sgp_avg.append(g2_sgp_v[i].sum()/len(g2_sgp))
        # calc stds.
        g1_std.append(g1_v[i].std())
        g2_std.append(g2_v[i].std())
        g1_sgp_std.append(g1_sgp_v[i].std())
        g2_sgp_std.append(g2_sgp_v[i].std())
    
    # average stds
    g1_avg_std = array(g1_std).sum()/len(g1_std)
    g2_avg_std = array(g2_std).sum()/len(g2_std)
    g1_sgp_avg_std = array(g1_sgp_std).sum()/len(g1_sgp_std)
    g2_sgp_avg_std = array(g2_sgp_std).sum()/len(g2_sgp_std)
        
    return [g1_avg_std, g2_avg_std, g1_sgp_avg_std, g2_sgp_avg_std]
            

def plot_conductances(rev_pot,std_avgs):
    subplot(len(pairs)/2+1,2,1)
    pylab.title("Conductances with various pairs\n of V_clamp recordings. (%3.1f)"%(rev_pot))
    for i in range(len(pairs)):
        j,k = pairs[i]
        subplot(len(pairs)/2+1,2,i+1)
        plot(t_imc,g1[i],t_imc,g2[i])
        imc_str = "_IMC"
        plot(t_sgp,g1_sgp[i],t_sgp,g2_sgp[i])
        sgp_str = "_SGP"
        pylab.ylabel("%2d & %2d"%(holding_V[j],holding_V[k]))
        pylab.legend(["ge"+imc_str,"gi"+imc_str,"ge"+sgp_str,"gi"+sgp_str])

    subplot(len(pairs)/2+1,2,len(pairs)+1)
    plot(t_imc,g1_avg,t_imc,g2_avg)
    plot(t_sgp,g1_sgp_avg,t_sgp,g2_sgp_avg)
    pylab.legend(["std = %4.2e"%(std_avgs[0]),\
                 "%4.2e"%(std_avgs[1]),\
                 "%4.2e"%(std_avgs[2]),\
                 "%4.2e"%(std_avgs[3])])

    show()

    

pairs = []
# we have (N-1)! pairs.
for i in range(len(holding_V)-1):
    for j in range(i+1,len(holding_V)):
        pairs.append((i,j))

g1     = range(len(pairs))
g2     = range(len(pairs))
g1_sgp = range(len(pairs))
g2_sgp = range(len(pairs))
g1_avg = []
g2_avg = []
g1_sgp_avg = []
g2_sgp_avg = []

std_avgs = calc_conductances(0.0)
plot_conductances(0.0,std_avgs)

